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ABSTRACT 

We present and describe IM3 SHAPE, a new publicly available galaxy shape measurement code 
for weak gravitational lensing shear. IM3 SHAPE performs a maximum likelihood fit of a bulge- 
plus-disc galaxy model to noisy images, incorporating an applied point spread function. We 
detail challenges faced and choices made in its design and implementation, and then discuss 
various limitations that affect this and other maximum likelihood methods. We assess the bias 
arising from fitting an incorrect galaxy model using simple noise-free images and find that it 
should not be a concern for current cosmic shear surveys. We test IM3 SHAPE on the GREAT08 
Challenge image simulations, and meet the requirements for upcoming cosmic shear surveys 
in the case that the simulations are encompassed by the fitted model, using a simple correction 
for image noise bias. For the fiducial branch of GREAT08 we obtain a negligible additive shear 
bias and sub-two percent level multiplicative bias, which is suitable for analysis of current 
surveys. We fall short of the sub-percent level requirement for upcoming surveys, which we 
attribute to a combination of noise bias and the mis-match between our galaxy model and 
the model used in the GREAT08 simulations. We meet the requirements for current surveys 
across all branches of GREAT08, except those with small or high noise galaxies, which we 
would cut from our analysis. Using the GREAT08 metric we we obtain a score of Q=717 for 
the usable branches, relative to the goal of Q=1000 for future experiments. The code is freely 
available from https://bitbucket.org/joezuntz/im3shape. 

Key words: methods: statistical, methods: data analysis, techniques: image processing, cos- 
mology: observations, gravitational lensing: weak, dark energy 



1 INTRODUCTION 

Weak gravitational lensing has the most potential for constraining 
the cosmological parameters that describe the dark universe (Al- 
brecht et al. 2006; Peacock & Schneider 2006). As dark energy 
slows the growth of dark matter structures, these in turn affect the 
appearance of distant galaxies, distorting them slightly as gravita- 
tional fields deflect their light. This distortion is known as cosmic- 
shear. Its measurement from individual galaxy images is crucial to 
the field of weak gravitational lensing, as the amount of distortion 
is directly related to the amount of intervening dark matter. How- 
ever, the stretch induced by cosmic shear is typically on the order 
of a few percent, which renders accurate shape measurement an 
intricate but also intriguing problem. 

Several upcoming and future observational surveys plan to 

* E-mail: jaz@star.ucl.ac.uk 



capitalise on the potential of cosmic shear for discovering the na- 
ture of dark energy. These include imminent ground-based projects 
(the Kilo-Degree Survey, Hyper Suprime-CatrQ and the Dark En- 
ergy Survey: DES3), ground-based telescopes under construction 
(the Large Synoptic Survey Telescope: LSST0), and future space 
telescopes {Eucli^ and WFIRST^). 

It is perhaps surprising that a major challenge for these exper- 
iments is simply measuring the shapes of the galaxy images. The 
difficulties arise because the galaxies involved are: small, less than 
the Point Spread Function (PSF) size in many cases; noisy, with a 
signal-to-noise ratio (S/N) of ~ 10; often irregular; and convolved 

1 |http: //www.nao j ■ org/Pro jects/HSC/HSCPro ject .html| 

2 http : //www . darkenergy survey . org 

3 http : / /www , Isst . org| 

4 http://sci.esa.int/euclid 

5 http://wfirst.gsfc.nasa.gov 
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with a PSF which itself induces ellipticity in the galaxy that is typ- 
ically greater than the gravitational shear we wish to measure. To 
fulfill the potential of cosmic shear these problems must be over- 
come, allowing precise and accurate measurements of the underly- 
ing gravitational shear. 

1.1 Existing methods 

The best known and most commonly-used m ethod in weak lens - 
ing shape estimation rema ins th at described by Kaise r et al. i cm, 
iLuppino & Raised J 19971) and iHoekstra et alj dl998l) . known as 
KSB or KSB+ (see also lSchneideJbOOd for an overview). KSB 
uses weighted quadrupole moments on galaxy images and PSF 
images to provide an estimator of gravitational shear. Many im- 
plementations of this approach, and derived methods adapted fo r 
space-based data (e.g. lRhodes et al.ll2000l . lSchrabback et alfcOld) . 
have been used to estimate gravitational shear over more than a 
decade. An important discovery of the Shear TEsting Programme 
(STEP) competitions was the sensitive dependence of the accuracy 
of K SB upon the details o f its implementation (see Section [L2l be- 
lowi lHevmans et alJl2006l : lMassev et alj2007h . 

Recent modifications to KSB have sought to improve ac- 
curacy by adding h igher-order terms to series approximations 
in the method (e.g lOkura & Futamasd [20091: Iviola et al.1 l201ll : 
iMelchior et all 1201 ll ; lOkura & Futamasell2012l) . More fundamen- 
tally different shape measurement approaches have also been pro- 
posed. Arelatively early example was the 'commutator' method of 
iKaised d200d) ; although this method has not found broad adoption 
in practice, the work contained a number of insights that were to 
inform the later development of the field. 

Gaussianization, or Re-Gaussianization, approaches 
(e.g., REGLENS: see iMandelbaum et alj I2012L Section 4.3; 



IMandelbaum et al I 120051 : lffirata& SeliaklHoolh first convolve 
galaxy images with an additional kernel designed to make the 
final PSF as Gaussian, and isotropic, as possible. Such images 
then closely match many of the implicit assumptions in weighted 
quadrupole moment approaches such as KSB, so moments can 
then be used to extract shear estimates with reduced bias. Gaus- 
sianization methods rema i n com petitive in terms of systematic bias 
( IMandelbaum et al1l2005l . l2012h and have been used successfully 
in a number of s tudies in th e Sloa n Digital Sky Survey (e.g . 



Huff et alj 



Hirata et al 



201 ll: iReves et all l20ld : IMandelbaum et alj 120061 : 



20041) . 



A family of methods known commonly as shapelets construct 
a model of PSF-convolved galaxy images as a sum of orthonor- 
mal Gauss-Lagurre polynomial functions (e.g. iBernstein & Jarvis 



2002MRefregier & Baconl2003l:lMassev & Refregieij2005[ iKuiikenl 
20061 : iMassev et alj 120071 : iNakaiima & Bernsteml koOTT) . Effec 



tively the ellipticity of derived shapelet mod els can then be used to 
provide an estimator of gravitational shear. IMelchior et al. I ( f20lch 
raised some conceptual concerns with the approach, and while sim- 
ilar (perhaps more drastic) conceptual problems with KSB did not 
affect its popularity there are fewer examples of sha pelet shear esti- 
matio n from real data in the literature (although see IVelander et al.1 

Honl). 

The number of proposed shape measurement methods is 
steadily growing. Interest in the problem has undoubtedly been 
stimulated by impending survey data of unprecedented statisti- 
cal power, and by the blind analysis contests we describe in 
Section [TT2l below. Examples of innovative recent proposals the 
competitions describe d in the next se ction, and include Fourier 
Domain Null Testing iBernsteirJ|2O10l) and shear estimation via 



lo okup table (e.g., 'Me gaLUT; see the review in the appendices 
of lKitching et alj|2012l) . 

The final family of shape measurement methods are com- 
monly known as model-fitting methods, and involve the compar- 
ison of simple parametric galaxy models (convolved by the PSF) to 
imaging data. The galaxy model ellipticity, among other parameters 
inferred in the mode l fit, make up t he she ar estimator. 

The lensfA l lMiller et al.1 120071 : iKitching et ail 120081 ; 
iMiller et"al] 1201 2h method seeks to perform Bayesian infer- 
ence on galaxy model parameters and performed well in the 
GRavitational lEnsing Acc uracy Teasting (GR EAT08) Challenge: 
see Section 11.21 below; dBridle et al.1 l2010l) . The more recent 
DeepZot method that also used a model-fitting algori thm proved 
the eventual winner of GREAT 1 Pitching et alj[2012l) . 

The im2shape public code dBridle et al.l2002lf! modelled both 
the ga laxy and the PSF as a sum of Gaussians, inspired bv lKuiikenl 
dl999h . It performed MCMC sampling and took the mean of the 
of the ellipti city values of the samples as the shear estimate. It 
was used by iBardeau et all d2005l 120071) to measure galaxy clus- 
ter masses. 



1.2 Shear measurement contests 

Researchers have come together to assess these many shear 
measurement methods in different regimes. The Shear TEst- 
ing Programme (STEP) was a joint blind analysis of simulated 
data by groups from within the shear m easurement community 
dHevmans et al]2006l : lMassev et al.l2007h . It showed that the shear 
measurement problem is far from trivial, but demonstrated that 
the cosmological constraints set at that time were not limited by 
shear measurement accuracy. The first two GRavitational lEns- 
ing Accuracy Testing (GREAT) Challenges posed the problem to 
computer scientists and provided sufficient simulatio ns to assess 
whether methods we re suitable for future surveys dBridle et al.1 
120081 : iKitchinelFioiol) . They showed that there are some regimes 
where existing methods work well enough, but others where fur - 
ther progress is necessary dBridle et al.l2010t IKitching et al.l2012h . 
The third challenge in the series, GREAT3Q is currently in progress 
towards a challenge opening and simulation data release in 2013. In 
addition to a large simulation dataset, GREAT3 is leading the col- 
laborative development of an open source extragalactic image sim- 
ulation software toolkit|f| further advancing the STEP and GREAT 
programme of providing the community with a common reference 
for comparison and improvement of galaxy shape measurement 
methods. 



1.3 Im3shape 

This paper presents a successor code to the im2shape code de- 
scribed above, called IM3SHAPE. It uses (by defaul t) a sum of 
two Sersic profiles for the galaxy, similar to lensfit dMiller et al.l 
|2007| ; IKitching et ai]|2008l ; IMiller et alJlIoTl) . and a Moffat pro- 
file for the PSF. It calculates the maximum likelihood parameter 
values and corrects for n o ise bias (see, e.g., Refregie r et alj|2012l; 
Melchior & Violal 1201 j Ifflrata et al.1 |2004 [Bernstein & Jarvii 




using the calibration scheme described in Kacprz ak et al.1 



6 |http : / /www . sarahbridle . net/ im2shape/| 

7 |http : /~7g reat3challenge . i nf o| 

8 https : / /git hub . com/ Gal Sim- developers /Gal Sim 
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Table 1. Shear measurement requirements on current, upcoming and future 
surveys. Requirements on the shear multiplicative and additive bias parame- 
ters nii and Ci, and the GREAT08 Q value are shown for survey parameters: 
area A, number of galaxies per square arcminute n ga i and median redshift 



Survey 


A 


"gal 


-111 


<», 


<H 


Q 


Current 


170 


12 


0.8 


0.02 


0.001 


43 


Upcoming 


5000 


12 


0.8 


0.004 


0.0006 


260 


Future 


2 x 10 4 


35 


0.9 


0.001 


0.0003 


990 



In Section [2] we discuss the requirements for shape measure- 
ment for current and future data sets. Then in Section[3]we present 
the IM3SHAPE code that we are releasing, the details of its opera- 
tion and the choice of parameters it uses. Section [4] is concerned 
with noise-free simulations, to verify the code in ideal circum- 
stances and to explore what biases the use of an incorrect model 
can introduce to the results. We present the results of tests of the 
code on the GREAT08 simulation set in Section[5] and conclude in 
Section[6] We give more details of the code implementation in Ap- 
pendix A, of the termination criteria in Appendix B and overview 
the analysis process in Appendix C. Appendix [D] gives the equa- 
tions and details used for analytic calculation of likelihood gradi- 
ents with respect to the parameters. 



3 IM3SHAPE METHODOLOGY & CODE 

The code we have developed and present in the remainder of this 
paper is called IM3SHAPE. It belongs to the family of model-fitting 
methods: it forward-fits a galaxy model to each data image it is 
supplied with and reports the parameters of the best fitting model, 
including the ellipticity components. For finding the best fit we use 
a maximum likelihood approach, i.e. we search for those model pa- 
rameters that minimize x 2 , the summed square difference between 
our generated model and the actual pixel data, weighted accord- 
ing to a user-supplied weight image. For details on the maxmimum 
likelihood estimation see Appendix lAl and |Pl 

Throughout we assume that we have selected a square portion 
of the image, with stamp_size pixels on a side. 

3.1 Code layout 

IM3SHAPE is a modular C code, with a significant amount of 
Python glue code that enables the setting up of new models and 
their options automatically. Key processes like optimization and 
galaxy model generation routines are encapsulated in complex 
structures with simple interfaces, rather than complicated func- 
tions. This architecture makes it particularly easy to design new 
galaxy or star models to be fit to data. Since the process of fitting to 
a large number of postage stamps is inherently trivial to parallelize, 
the core of the code is all single-threaded. 



2 REQUIREMENTS 

Following iHevmans etal] l2006t) . the accuracy of shear measure- 
ment methods are generally quantified in the literature in terms of 
the multiplicative error m and additive error c on the true shear 7*, 
such that 



7i = (1 + m,i)jl + a 



(1) 



where 7, is the shear estimate and the subscript i refers to the two 
shear components. It is assumed there is no cross contamination 
between, for example, 71 and 72. 

The requirements on m a nd c for cosmic shear are computed 
in Amara &R6fregieil d2008l) as a function of the survey area, 
depth and galaxy number density (see also iKitching et aill2009r) . 
lAmara & Refregierl d2008[) consider a two-point statistical analysis 
of the shear field, and the requirements are set so that the system- 
atic error equals the statistical error, thus providing an upper limit 
on the allowed bias. 

The GREAT08 feridle et al.l2008h challenge quantifies the ac- 
curacy of shear measurement methods using a single number called 
th e quality factor, Q. Th e equation relating Q to m and c is given 
in IVoigt & Bridie] dioioh. A related but different definition is used 
in GREAT 10 dKitching eUu"ll2012ri but we stick to the GREAT08 
definition here because we use the GREAT08 simulations. 

In Table Q] we show the requirements on m, c and Q 
for three sets of survey parameters, chosen to represent cur- 
rent, upcoming and future cosmic shear surveys. Examples 
of surveys with requirements comparable to the 'current' set 
of requirements are the Canada-France-Hawaii Telescope Lens- 
ing Survey (CFHTLenS: IHevmans et al.l 1 20121), and likely the 
first results from KIDS, PES and HSC jde Jong et alj 120121 : 



iThe Dark Energy Survey Collaboration! 120051: lTakadall2010l) . The 
upcoming survey requirements correspond roughly to the full DES 
and HSC surveys. Th e future survey requirements correspon d to 
LSST and Euclid (e.g. lChang et alJl2012l : lLaureiis et al.11201 ll) . 



3.2 Model choice 

Using the flexible systems of IM3SHAPE we have experimented 
with fitting a variety of models, mainly to the images in the 
GREAT08 challenge. 

We chose by default a two-component Sersic, with a de Vau- 
couleurs bulge and an exponential disc, with the same centroid and 
ellipticities. In addition we simplify the fitting by fixing the ratio of 
scale radii, at unity[f| 

Within that general model we can choose different parame- 
terizations for some components - these are discussed in detail in 
Appendix [A] 

Our complete set of parameters is shown in Table [2] The pa- 
rameters marked as fixed were kept at default values during our 
standard analysis but can be varied easily if desired or to study 
model bias. 



3.3 Model fitting 

Generating the model image described in the previous Section is 
not trivial: we need to ensure that the model is at high enough reso- 
lution to capture sharp central galaxy features. This is discussed in 
Section|4~T1 

Assuming images of the galaxy model can be accurately gen- 
erated (see Section |4~H , maximizing the computed likelihood (i.e. 
minimizing the weighted, squared difference between the non- 
linear model and the data) is technically and computationally chal- 
lenging. Reaching the current level of speed and accuracy in max- 
imization within IM3SHAPE required considerable exploration of 



9 This is similar but not identical to the model adopted in the GREAT08 
simulations. IM3 SHAPE should therefore work reasonably well on those im- 
ages, and we note that this model was originally chosen to be a reasonable 
description of real galaxies. 
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Parameter Meaning Fixed 





Horizontal centroid 




yo 


Vertical centroid 




ei 


x-y shear 




'■2 


45° shear 




rd 


Disc half-light radius 




A b 


Bulge peak flux 




A,, 


Disc peak flux 




R r 


Bulge-disc size ratio 


/ 


rid 


Disc Sersic index 


/ 


m, 


Bulge Sersic index 


/ 


Ae 


Bulge-disc ellipticity 


/ 


AO 


Bulge-disc angle 


/ 



Table 2. The complete set of parameters used in the default galaxy model. 
Fixed parameters were not by default varied during parameter optimization. 
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Figure 1. Slices in minus log-likelihood in six of the principal parameters 
(see section [A} around the best fit point found for one galaxy by the likeli- 
hood maximizer, which is shown by a star. 



maximization algorithms and their settings. We therefore give a 
brief overview of our findings in Appendix [B] and an overview of 
the whole process of image generation and fitting in Appendix Icl 

With tuned parameter settings the maximizer, which uses a 
Levenberg-Marquardt algorithm, works well. Figure Q] shows ex- 
ample likelihoods from the first galaxy of set 3 in GREAT08. It 
can be seen that the point is indeed a good minimum (and further 
searches can confirm it to be a global one). The asymmetry in some 
of the parameters is pa rt of the cause of the noise bias phenomenon 
dKacprzak etail 20121) . 



4 NOISE-FREE SIMULATIONS 

Generating the model image described in the previous Section is 
not trivial; for example, we need to ensure that the model is at high 



enough resolution to capture sharp central galaxy features without 
aliasing. We also need to understand how best-fitting models may 
be biased in cases where they cannot represent the data perfectly. 
In this Section we investigate the behaviour and performance of 
IM3SHAPE on noise-free simulations to find optimal parameter set- 
tings for model generation, and analyse in a simple way what po- 
tential biases can arise when the model is insufficiently flexible to 
describe the data. 

4.1 Upsampling biases 

The IM3SHAPE software uses the Discrete Fourier Transform 
(DFT) to render images of convolved galaxy profiles: this can 
be done with speed thanks to the Fast Fourier Transform (FFT) 
algorithm, freely avail able in highly efficient implementations 
dFrigo & Johnsonll2005l) . However, representing the convolution of 
telescope PSFs with analytic, continuous, non-bandlimited func- 
tions such as the Sersic profile requires careful numerical approxi- 
mation, as can (if the PSF is also not bandlimited) the subsequent 
rendering of convolved profiles into pixellated images (see, e.g., 
lMarksll2009l) . Decisions must be made about how and where to 
make approximations in the representation of these profiles, bal- 
ancing precision against computational cost. 

Upsampling, which is the more accurate rendering of model 
images using higher resolution (upsampled) intermed iate images 
for in ternal operations such as convolution (see, e.g., iRowe et al .1 
is a key tool to eliminate defects in modelling due to alias- 
ing. We therefore investigate the requirements on the upsampling 
parameters used by the image generation code in IM3SHAPE. These 
three resolution parameters are: 

( 1 ) Ui , called upsampling in parameter files - the convolution be- 
tween the galaxy model and the PSF model is performed on a 
grid of (Ui x stamp .size) 2 pixels in size, i.e. each image pixel 
is divided into l]\ sub-pixels across. The flux within each sub- 
pixel is computed assuming the intensity is constant within it 
and equal to the intensity at the centre of the sub-pixel, for all 
pixels except the central N 2 pixels (see (2)). 

(2) N u , called njpixels Jo jup sample in parameter files - within a 
central galaxy region of N 2 image sub-pixels the flux is inte- 
grated on a even finer grid within each sub-pixel (instead of 
assuming a constant intensity across each sub-pixel). 

(3) U2, called n.central .pixel jupsampling in parameter files 
- determines the number of integration intervals within a sub- 
pixel in the central image region; therefore the overall resolu- 
tion in the central galaxy region is Ui x U2 times the data image 
resolution. 

The required values of the these resolution parameters depend on 
the survey (see section |2j- In figure [2] we plot the multiplicative 
bias as a function of the Ui parameter for different values of N u 
and with U2 = 7. (The additive bias is negligible in comparison 
to the multiplicative bias because we use a circular PSF for this 
investigation.) Results are shown for a pure bulge and a pure disc 
with Rgp/Rp — 1.4 and 1.14, which are sim ilar to the fiducial 
and small size branches used in GREAT08 (see lBridle et al.|[2010t 
incl uding for size definition). The bias is computed using a ring- 
test dNakaiima & Bernsteinll2007l) with 4 angles, {0, 45, 90, 135} 
degrees. 

The multiplicative bias is obtained by shearing each galaxy in 
the ring by an amount 71 = 0.02 along the a>axis and 72 = 0.01 
along the line y = x. The PSF is a circular Moffat profile with 
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Figure 2. Resolution requirements: multiplicative bias (mi crosses, ni2 open circles) as a function of the number of sub-pixels U\ used in the convolution 
between the galaxy model and the PSF model. Black solid, red dashed, yellow dot-dashed and blue dotted lines are for integration of the galaxy flux within 7x7, 
5x5, 3x3 and lxl sub-pixels around the centre of the galaxy. The flux is integrated over a 7x7 grid within each sub-pixel. Top panels: small (R gp /R p = 1.14) 
galaxies. Bottom panels: fiducial size (R gp /R p = 1.4) galaxies. Left panels: de Vaucouleurs galaxies. Right panels: exponential galaxies. The shaded regions 
show the requirements on rrii given in Table [T] The upper edge of each shaded region (from bottom to top) shows the upper limit on the bias allowed for 
future, upcoming and current surveys respectively. 



FWMH=2.85 pixels and shape parameter /3 = 3.5. The postage 
stamp size (stampsize in parameter files) is 25 pixels. We find that 
the requirements on the resolution parameters do not change if we 
increase the stamp size to 45 pixels. The bulge-to-disk scale size 
ratio, R r , is fixed to 1 during the fitting. 

We find that for an upcoming survey (medium green shaded 
area in figure [2} we require Ui J?5, N u ^ 3, and U2 J?7. For 
U% < 7 the results are unstable. 



only touch upon these iss ues here, we refer the i nterested reader to 
IVoigt & Bridie] feOld) and lMelchior et al.l d2010l) for a more exten- 
sive disucssion. 

Except where otherwise noted, our galaxy simulations in this 
section used R gp /R p — 1.4, equal component scale radii, and 
half the flux in each of the bulge and disc components Fb / {Fb + 
F D ) = 0.5. The PSF is a Moffat with ei = e 2 = 0.03, ,8 = 3, and 
FWHM= 2.85. 



4.2 Model bias 

Having established that the resolution parameters described above 
are high enough to remove any bias from aliasing effects, we now 
turn to the issue of model or external bias arising from fitting an 
incorrect model to the data. If a forward model approach could 
be petfectly implemented and the posterior distribution accurately 
characterized for each galaxy (rather than being reduced to a single 
point-estimate like the maximum likelihood) and combined appro- 
priately, then this would be the only remaining source of bias for 
simple images like those in GREAT08. In practice, there may be 
other problems such as neighboring galaxies and image artefacts. 

Very realistic galaxy models are almost certainly out of reach 
of forward methods for the forseeable future - modelling spiral 
arms, for example, is hard enough for a single high-resolution im- 
age, let alone 10 8 noisy ones. We defer testing on extremely re- 
alistic models (for exam ple using the SHERA simulation code of 
Mandelbaum et al .120121 . or Galsim), and instead apply a necessary 
but not sufficient test - fitting a simple but incorrect model using 
IM3SHAPE. 

Again, we perform these tests without noise, to separate the is- 
sue of noise bias on maximum-likelihood estimates from external 
bias; of course, there may be an interaction between the two. As we 



4.2.1 Sersic Indices 

We first test a simple difference between the input (true) model and 
the fitted one - a different index of the Sersic profile. Our input 
(true) model is generated with a single component, a Sersic profile 
with index n s , however, we fit our standard bulge plus disc model. 

Figure[3]demonstrates this bias: it shows the pre- and post-PSF 
convolved simulated galaxy images, and the residuals between the 
simulated galaxy image and best-fit galaxy model, for two different 
values of the true Sersic index. The images are noise-free and the 
correct PSF model has been used in the fits. The error on the true 
galaxy ellipticity is approximately 10~ 3 for n B — 1.7, and 2x 10~ 4 
for n 3 ~ 3.5. 

The residual image shows different behaviour along the x 
and y directions, which is what causes the ellip ticity to be mis- 
estimated. As discussed in lVoigt & Bridle! d2010l) . the galaxy size 
is worst estimated along the direction where the galaxy is smallest 
relative to the PSF. Therefore the ellipticity is incorrect. 

Our model used in IM3SHAPE is correct when n 3 = 1 or 
n 3 = 4, since at these points the true model is a subset of the 
IM3SHAPE model. We generate simulations with our fiducial sets 
of galaxy parameters, except for the Sersic index, which we vary. 
We apply a ring test to the simulations to assess the multiplicative 
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Figure 3. Examples of model bias: IM3SHAPE co-elliptical bulge (n s = 4.0) plus disc (n s = 1.0) model fits to simulated single-component Sersic galaxies 
with n s = 1.7 (top panels) and n s = 3.5 (bottom panels) . Postage stamp images are shown for the simulated galaxy (left-hand panels), the PSF-convolved 
simulated galaxy (middle panels), and the residuals between the simulated galaxy and the best-fit galaxy model from IM3SHAPE (right-hand panels). The 
simulated galaxy has ellipticity components ei = 0.2 and e2 = 0.0 and the radius along the major axis is 1.84 pixels. The PSF is a circular Moffat profile 
with p = 3.0 and FWHM=2.85 pixels. 




True Sersic Index 



Figure 4. The mean multiplicative bias m as a function of true model Sersic 
index when fit with a bulge+disc model, for three different image sizes as 
shown. The shaded bands are the current and upcoming survey target biases. 
Since the simulations were noise-free error bars are negligible. 



m and additive c biases; the former is shown in figure |4] a figure 
for the latter is omitted as it is negligible in the noise-free case. 
As expected, the bias becomes zero when the true and fitted model 
coincide at n s — 1 and n 3 — 4. 

Figure[4]suggests that this kind of model bias has the potential 
to be a problem for upcoming surveys but is only marginally rele- 
vant for current generation of experiments. Allowing for additional 
flexibility of the model could resolve this potential issue. Further 
tests are needed on more realistic galaxy images to assess the true 
impact on upcoming surveys. 

4.2.2 Different bulge and disc ellipticity 

In figure [5] the multiplicative bias arising from a slightly more re- 
alistic model difference is shown - a bulge and disc model with 
components that are not co-elliptical. This is of interest because it 



is a difference between the IM3SHAPE model and the GREAT08 
simulations. The input image bulge ellipticitsP'l e& is fixed at 0.3 
and the disc ellipticity varied so that et, — + Ae. The difference 
between the bulge and disc ellipticity, Ae, is varied between -0.25 
and 0.25. Typically, bulges are rounder than discs, i.e Ae < 0. The 
model fit to the data is our standard two-component co-elliptical 
profile. 

Once again a ring test is applied to assess the bias, and the 
noise is zero. The other galaxy parameters are set to the fiducial 
values described above, except that the different curves show dif- 
ferent relative amplitudes of the bulge and disc components. The 
bias is zero when either Ae is small or the galaxy is completely 
bulge- or disc-dominated, since in that case the secondary compo- 
nent does not affect the fit at all. The effect is greatest when the 
components have comparable total fluxes, and this bias too is small 
for the current generation of experiments, if bulges are rounder than 
discs. 



4.2.3 Incorrect Radius Ratio 

As depicted in Table(2] the ratio between the half-light radii of the 
bulge and the disc is fixed by default in the code. This is partly 
to reduce the number of parameters (since fewer parameters gen- 
erally means easier optimization), but also because this particular 
parameter causes problems with convergence in some regions of its 
range. 

This is clearly a potential source of model bias - the ratio of 
these radii varies in real galaxies (when bulges and discs can be 
clearly identified). Furthermore, in GREAT08 this value is around 
R r ~ 1.5. We again apply a ring test to noise-free simulations 
to check how severe this problem is. If we alter the ratio without 
changing the disc radius, then the bulge radius changes accordingly. 



10 We define ellipticity as e 
Appendix lAl 
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Figure 5. The mean multiplicative bias m as a function of true model 
Ae, the difference between bulge and disc ellipticities, when fit with a co- 
elliptical bulge+disc model. The shaded bands are target biases for current 
and upcoming experiments. Since the simulations were noise-free error bars 
are negligible. In real galaxies Ae is usually negative. 
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Figure 6. The mean multiplicative bias m as a function of the true ratio be- 
tween the bulge and disc half-light radii for various galaxy sizes is shown. 
Note that the bias arises as the ratio has been fixed to one during the fitting 
with IM3SHAPE. Again, the simulations are noise-free so there are no error 
bars. The shaded bands are target biases for current and upcoming experi- 
ments. 



To separate the effect of fitting small galaxies from this model bias 
we keep the ratio of the convolved image to PSF size, i.e. R gp / R p , 
fixed at three different values corresponding to GREAT08 sizes. 

Figure [6] shows the results of this test. If real galaxies have 
smaller bulge scale radii than discs then the biases are mostly 
smaller than the requirements for upcoming surveys. If they are 
larger then the bias would be acceptable for current surveys but 
problematic for upcoming surveys, with a sign which depends on 
the galaxy size. For the GREAT08 fiducial branch (R — 1.4, true 
R r ~ 1.5) we expect a slight over-estimate of the shear. This will 
combine with the slight underestimate expected from figure 5 for 
bulges which are rounder than discs, and may partially cancel out to 
bring results down to be compatible with the requirements for up- 
coming surveys. We recommend a program of testing on realistic 
images for calibration of upcoming surveys. 



5 GREAT08 

We tested IM3SHAPE on the Gravi tational LEnsing Accurac y Test- 
ing 2008 (GREAT08) Challenge dBridle et al.|[2008ll2010h which 
has a sufficient number of galaxies to test methods at the accuracy 
required for future surveys (unlike STEP1 and STEP2). Further- 
more, it features a constant shear across the images which greatly 
facilitates the computation of success metrics (unlike GREAT10) 
and therfore performance evaluation. The GREAT3 Challenge is 
still under construction (Mandelbaum, Rowe et al. in prep). 

The GREAT08 Challenge includes two blind competitions: 
one with low noise, which has a signal-to-noise ratiq 11 ! of 200 
('LowNoise_Blind') and one with more realistic signal-to-noise ra- 
tio around 20 ('RealNoise_Blind'). In this section we show the re- 
sults of IM3SHAPE on these simulations. 



5.1 Low Noise Blind 

The LowNoise_Blind challenge consists of 15 image files, each of 
which has 10,000 galaxy postage stamps of 39x39 pixels. All the 
galaxies are a mixture of bulge and disk components which are co- 
centered but not co-aligned. The misalignment angle is drawn from 
a Gaussian with a standard deviation of 20 degrees. We therefore 
expect a model bias to apply in these cases as described in sec- 
tion |4.2| The sig nal-to-noise ratio of 200 is large enough that noise 
bias described in Ka cprzak et al .1 l l2012h is expected to be negligi- 
bl e. The PSF and ga laxy size parameters are described in detail 
in Bri dle et al. lima. In summary, there are 5 image files of each 
of 3 different galaxy sizes. This allows us to show how the success 
metrics depend on galaxy size. 

The GREAT08 images are simulated at essentially infinite res- 
olution by using p hoton shooting for the PSF and galaxy profiles 
dBridle et al]|2010h . IM3SHAPE makes its models at high resolu- 
tion (see section |4~l"l l. For this test we ran at two different values 
of this resolution, with an upsampling Ui = 5 for a coarse run 
and Ui = 7 for an improved, higher resolution run that took about 
twice as long to complete. We are interested to see both results on 
the small LowNoise_Blind dataset and would like to use the coarser 
resolution for the much larger RealNoise .Blind dataset. 

In figure [7] we show the IM3SHAPE results on the 
LowNoise_Blind GREAT08 data. At high resolution the resulting 
Q-factor scores reach level Q > 1000 for all branches, a level at 
which the statistical noise in the challenge dominates - scores above 
this level are consistent with zero error. We reach this score in spite 
of a model bias: the the GREAT08 simulated galaxies have mis- 
aligned isophotes, whereas we fit co-elliptical ones, and the ratio of 
bulge to disc size ratio is different. 

At the coarser resolution we are within the requirements for 
future surveys for all but the multiplicative bias on small galax- 
ies, which is nonetheless acceptable for the current generation. We 
therefore use the coarser resolution in the remainder of this paper. 

We also show the results submitted by others at the time 
of GREAT08. Significant developments have been made in other 
codes in the meantime, but we can conclude that at that time there 
were no methods as effective as IM3SHAPE in this high signal-to- 
noise regime. The closest competitor is the stacking method by AL, 
which has not yet been developed into a method that can be used 
on realistic data. We should however note that the model in the 



11 Here, we ado pted the GREAT08 definition of signal-to-noise ratio de- 
fined by (Al 1) in lBridle et all l2O10l) , Appendix A, page 14. 
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IM3SHAPE code is very well matched to the GREAT08 simula- 
tions, in that they both use exponential disc plus de Vaucouleurs 
bulge galaxies. An investigation into this topic is beyond the scope 
of this work, and will be addressed by future challenges including 
GREAT3. 

In running this analysis we determined which parameters in 
the code can be tweaked to minimize runtime while retaining high 
accuracy. The most important such parameters are those controlling 
the details of the LevMar optimization and convergence behavior. 
In subsequent work we have found that tuning these parameters is 
quite specific to the noise levels and other parameters of the galax- 
ies in question, and we have not been able to draw any general 
conclusions about them. 



5.2 Real Noise Blind 

The GREAT08 RealNoise_Blind Challenge is considerably larger, 
with 2700 sets of 100 x 100 galaxy stamps. The 2700 sets are di- 
vided into 9 branches which span a range of simple galaxy mod- 
els and PSF ellipticities and a range of realistic galaxy sizes and 
noise levels. The fiducial branch has a Signal-to-Noise a factor ten 
lower than in LowNoise_Blind but the same galaxy size as the cen- 
tral branch of LowNoise_Blind and the same PSF and galaxy type 
(co-centred but not co-elliptical bulge plus disc). Two branches ex- 
plore the same two extrema in galaxy size as LowNoise_Blind; two 
branches vary the Signal-to-Noise by a factor either side of the fidu- 
cial value; two branches vary the PSF ellipticity by a factor of two 
and in direction; and two branches explore variants on the galaxy 
model: bulge or disc and bulge plus di sc with an offset in t he center. 
For more details on the challenge, see lBridle et al.ld2008h . 

We tuned the maximization parameters described in Appendix 
IB II to improve the mean likelihood of all the galaxies in a single 
set; on the basis that if a set of minimisation parameters did not 
find the best likelihood then it must not have converged fully. We 
note that the mean ellipticity obtained for this single set was quite 
sensitive to the maximization parameters, and therefore we recom- 
mend users of the IM3SHAPE code should perform a similar tuning 
on new galaxy samples. 



5.3 Noise Bias Calibration 

At the noise levels in the RealNoise_Blind Challenge we do ex- 
pe ct that we need to co rr ect for the nois e bias effect described 
in iRefregier et al.l d2012h . IKacprzak et al.l J2012h and references 
therein. We describe our approach to noise bias calibration in this 
subsect ion. 

In IKacprzak et alj l2012h we ran IM3SHAPE on noisy simu- 
lations with a range of input parameters to find the behaviour of 
the noise bias with galaxy properties. We then attempted to remove 
the bias by matching the measured galaxy properties to the simu- 
lation input parameters to predict the expected bias, which we sub- 
tracted off. We ran this method on branch 3 of our GREAT08 anal- 
ysis, which has the same galaxy model as IM3SHAPE, and found 
Q = 166, with m = 3.0 x 10~ 2 and c = 7 x 10" 4 which is 
nearly satisfactory for current experiments but well outside our re- 
quirements for upcoming experiments. 

This straightforward correction fails to do better because the 
observed properties are themselves noisy (and thus noise-biased), 
and so our bias estimate is itself biased. In practice, we have 
found it ineffective to fully calibrate galaxy ellipticity measure- 
ments one-by-one using their individual properties, and turn instead 



to a population-based approach in which a complete class of galax- 
ies is calibrated collectively. 

The approach that we take is discussed in Appendix|E] Briefly, 
we use the deep data from the low-noise part of the challenge to 
provide distributions of the galaxy properties, and then perform 
simulations to compute the mean multiplicative and additive biases 
for those properties. We then apply these mean biases to the shear 
estimates from each set. 

To predict biases for a given parameter set we simulate galax- 
ies at grid points with scale radii Vd = {1,2.5,5} and bulge flux 
fractions F = {0, 0.5, 1}, and with ellipticity from ei = —0.8 
to 0.8. At each grid point we fit a cubic polynomial to the bias 
e — etruc as a function of e tr uc, and store the coefficients. To get 
the bias for intermediate parameters we use a Gaussian radial basis 
function interpolator. Since the simulations are costly we simulate 
only at S/N = 20, and ext rapo late to other values usi ng the result 
from lRefregier et al1d2012l) and lKacprzak et al!d2012l) that the bias 
scales as the inverse of the square of the signal-to-noise. 

There are several factors about the structure of GREAT08 that 
are unrealistic in this context. Firstly, not every branch in the main 
sample has a corresponding low-noise version, whereas in real sur- 
veys it should be possible to match populations with more care in 
most cases. This would tend to make noise calibration harder on 
GREAT08. Secondly, the sizes of galaxies in GREAT08 are not 
drawn from a population - they have a single true value. This would 
tend to make calibration more effective than is realistic. We also 
note that the information needed to fully perform the calibration in 
this way was not public at the time of the challenge. 

Where possible we match branches using the corresponding 
low-noise branch; where none exists we use the nearest approxi- 
mation. See Appendix lElfor more details. 



5.4 Scores 

Q-factor scores and m and c values for the pre- and post-noise bias 
calibration results are shown in figure [8] Even before the calibra- 
tion the bias on the S/N = 40 galaxies was small enough that we 
achieved a quality factor Q > 1000, with m and c small enough 
for upcoming surveys. At lower signal to noise, as we expect, the 
pre-calibration Q factor drops to below 100. 

After the noise bias calibration procedure described above the 
Q factors and the c values reach the levels required for upcoming 
experiments for most of the branches of the challenge. The results 
are stable with PSF type. The m values are somewhat more variable 
but, depending on the exact branch, can reach the required levels. 
In particular the value is good in the b/d branch, where where we 
expect no model bias. 

Our most encouraging result is the stability across the galaxy 
type branches. The only branch for which we use the correct model 
is the first galaxy type, where the simulation model is single- 
component bulge or disc. For all other branches we expect some 
model somewhat like the ones discussed in section |4~2l but as in 
that section the effect is not critical for current surveys. 

The branches for small and noisy galaxies are clearly more 
problematic, and we would be forced to remove such galaxies from 
any current analysis. We believe the extrapolation fails badly in the 
high-noise case because of a slowdown in the magnitude of the 
bias with ellipticity, which does not follow the b ~ (S/N)~ 2 re- 
lation that we used. Very noisy galaxies have much broader scatter 
in measured ellipticity, and when these values start to push up to 
the edge of the space, at |e| = 1, this acts as an m < bias which 
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Figure 8. The Quality factor score Q, multiplicative bias m, and additive bias c obtained in the RealNoise-Blind section of the Great()8 challenge. Each 
column shows scores with a different variation from the fiducial galaxy parameters. The thick solid blue line is our primary result, the IM3SHAPE scores 
following the matched-population noise bias calibration described in the text. The dashed green line is the score without any calibration. The green bands are 
the requirements for current and upcoming surveys. The thin grey lines are other entrants at the time of the challenge, with dots indicating stacking methods. 
The galaxy types are bulge or disc (50% of each; labelled b/d), co-elliptical bulge plus disc (b + d) and bulge plus disc with offset centroids (oft). The PSF 
values were fiducial (fid), rotated by 90°, and double ellipticity(e X 2). 



partially counteracts the usual m > one. Missing this effect can 
lead to extrapolated values over-correcting the bias. 

6 CONCLUSIONS 

We have presented a public galaxy shape fitting code, IM3SHAPE, 
for weak lensing shear measurement. Our code efficiently finds the 



maximum likelihood fit of our model to any given galaxy, and does 
so extremely well: we find little sign of any bias due to poor min- 
imization for simple simulated galaxies, following careful tuning 
of the model image resolution settings and other parameters. The 
code is modular and easily extensible, and we would welcome com- 
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munity contribution to any part of i The technical problems in 
generating models and fitting them to galaxies are extensive, with 
particular problems of resolution, parameterization, and non-linear 
numerical optimization. We discuss our solutions to these problems 
in the appendices to this paper. 

We find that we need to generate model images at five times 
the data resolution to reach our desired accuracy for upcoming sur- 
veys. We also find that to accurately capture small scale features in 
de Vaucouleurs cores we must sample a small central region even 
more finely - at 35 times the data resolution. 

Our code fits a simple bulge plus disc galaxy model to galaxy 
images, and we have made an exploration of model bias, which can 
arise when true galaxies fail to be as neatly analytic in form as our 
model ones. We find that reasonable deviations between the sim- 
ple true and assumed models cause only minor bias, not significant 
for current data and early results from the upcoming generation of 
surveys. The bias will, however, attain greater significance as up- 
coming surveys reach their full statistical potential. Our tests were 
rather limited in that we considered only simple differences from 
our canonical model - deviations in Sersic index, and differences 
between the size and ellipticity of the two model components. We 
look forward to performing more realistic tests with the Gal Sim 
tool set in future work (Mandelbaum, Rowe et al. in prep). 

We have demonstrated the utility of our method and code 
by testing on the GRaviational lEnsing Accuracy Testing 2008 
(GREAT08) Challenge in both low and realistic noise regimes. We 
find extremely good results for the low noise simulations, with 
multiplicative and additive shear measurement biases well within 
the requirements for upcoming (Stage III) surveys, for the medium 
sized and large galaxies. For practical reasons we chose to compro- 
mise on the resolution we used in the fits and thus our results for 
the sma llest galaxies are only s ufficiently good for current surveys. 

In iKacprzak et alj j2012h we calculated the degradation of 
shear measurements due to noise on images and discussed correc- 
tion schemes. We do not expect this to be significant at the signal- 
to-noise of the LowNoise simulations (S/N = 200). We expect the 
main difficulty to arise from resolution and the differences between 
our fitted model and the GREAT08 galaxy models. Both our galaxy 
models and the GREAT08 low noise simulations assume galaxies 
are made of a de Vaucouleurs bulge plus an exponential disc. The 
difference is that we assume these two components are co-centric 
and co-elliptical with the same size ratio, whereas in GREAT08 
low noise the two components have different ellipticities, are not 
co-aligned and have a bulge scale radius which is approximately 
50% larger than the bulge scale radius. We find that there is little 
resulting problem arising from model bias alone, as expected from 
our simpler investigations. 

We ran IM3SHAPE on the 27 million galaxies in the realis- 
tic noise part of the GREAT08 Challenge, which we chose over 
GREAT 10 because it is simpler and the variable shear fields used 
in the latter were not relevant to this method. 

At this signal-to-noise (S/N ~ 20) w e do expect significan t 
noise bias due to pixel noise in the images dKacprzak et alj|2~012h , 
leading to a multiplicative shear measurement err or of a few per 
cent. We found that the correction scheme tested in dKacprzak et al .1 
120121) does not work well on GREAT08 due to the bias-on-bias 
problem described in that work. We therefore deviated from the 
methods discussed there, which corrected each galaxy individu- 
ally, in favour of a global correction using the knowledge of galaxy 



properties measured in the corresponding low noise sample. The 
leading order correction to shear power spectrum measurements is 
expected to simply depend on the mean multiplicative bias of a 
population, which can be corrected globally. Since the noise bias 
depends on galaxy properties like size, and galaxies with similar 
properties can be clustered, then it is possible this scheme may not 
be satisfactory for high precision corrections. 

We scale the noise bias correction from the fiducial value ac- 
co rding to signal- t o-nois e, as expected in the analytic calculations 
of iRefregier et all d2012l) . This works well for the lower noise im- 
ag es but less well for th ose at higher noise levels as also noted 
in dKacprzak etai]|2012l) . As stated in that work we are not sat- 
isfied with the performance of our optimiser in this highest noise 
regime and would cut such images from cosmic shear analyses 
at this point, or calibrate using a more extensive suite of popula- 
tion matched simulations. A similar concern applies to the smallest 
galaxies. 

For all other branches of GREAT08 we attain additive and 
multiplicative biases which are within the requirements for current 
surveys, and have quality factor metrics which are as good or bet- 
ter than other non-stacking methods at the time of the GREAT08 
Challenge. 

Further work appears to be required before IM3SHAPE can be 
applied to the final years of upcoming Stage III cosmic shear sur- 
veys. We speculate that our residual shear measurement bias is due 
to the combined effects of using the incorrect galaxy model in the 
fit, and the noise on the images. Ideally an improved shear measure- 
ment method could be found that did not suffer from these types of 
bias. Alternatively the calibration scheme needs to be tailored more 
accuractly to the analysis in question, for example by using a suite 
of low noise images with known shears. 

To be fully applicable to real data there are a number of ex- 
tensions required to the code, for example (i) the ability to handle 
multiple exposures of the same patch of sky (ii) the capacity to deal 
with or flag neighboring objects contaminating the shear measure- 
ment. 

Furthermore, we have not investigated here the effect of 
postage stamp size or the effect of a weight map cutting out parts 
of the galaxy image. These developments are ongoing but testing 
them is beyond the scope of the present work. 
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APPENDIX A: MODELLING AND PARAMETERIZATION 
Al Choice of ellipticity parameters 

The performance of all the optimizers we have tried is heavily de- 
pendent on the particular choice of parameterization. In particular, 
the choice of how to parameterize the shear components is partic- 
ularly important since they are our target. We define ellipticity to 
match the effect of shear on images: e = (a — b)/(a + b). The two 
obvious parameterizations based on this, (ei, e2) and (e, 8), were 
both found to be problematic. 

Using the two separate ellipticities ei and e2 as parameters 
creates a non-linear boundary at e? + e| = 1, the edge of the unit 
circle. 

While most optimization methods cleanly handle linear 
boundaries (for example simple limits on parameters), non- linear 
ones are less well served. We tried simply assigning an extremely 
low likelihood to any proposed sample beyond this boundary, but 
fast optimization methods use likelihood differences to determine 
step sizes, so they typically behave wildly if large negative values 
are used to signal edges. Since the images of model galaxies di- 
verge as |e| — > 1 simply making the models of these galaxies is 
also prone to numerical problems. 

Using the polar coordinates (e, 9) as parameters works well at 
high ellipticity where we use a simple parameter limit on e, but fails 
at lower ellipticities. As e — s> the likelihood becomes a plateau 
in 6 since the coordinates are degenerate. Optimizers spent an ex- 
tremely long time wandering around the circle. 

We found a reliable parameter set which is a variant of the 
first parameterisation, (ei, e^). We use (ei, e-i) as the parameters 
varied by the optimizer, but map them to a disc wherever they stray 
beyond it. We mirror the parameter space radially about some e 
just below unity (we found 0.95 worked well), so that as the input e 
increases above the maximum, the value given to the minimiser first 
decreases and eventually reaches zero, and turns negative so that 
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the signs of the output e\ and e2 are flipped. Although this mirror 
structure causes repetition in the likelihood surface and multiple 
modes we have not found this to be an issue in practice. 

A 1.1 Other parameters 

There are various radius parameters we could choose for each of 
the components. The choice is particularly important for the sharp 
bulge profiles, as their full-width half-maxima are extremely small 
and so they are not suitable parameters. 

A better choice that works for both components is the radius 
which encloses half the total flux of the complete profile. Rather 
than using two radius parameters, we the disc radius and the ratio 
between the radius of the bulge and disc. We find (see section |4~2l 
that no large model bias is generated by fixing the latter parame- 
ter provided it is approximately correct. We therefore use a fixed 
value for any given data set. For our GREAT08 runs this value was 
arbitrarily set to unity, to avoid trying to match too closely to the 
GREAT08 truth values. 

A likelihood plateau occurs for small objects like stars, where 
the shears are irrelevant. In these regimes the fitted radii of the 
objects tended to fail to converge as they varied around zero. We 
solved this problem by setting a small transition radius R — 0.4 
sub-pixels and a replacing the radius r with \(R + r 2 / R) when 
\r\ < R and \r\ otherwise. The minimum value radius actually 
used is then R/2 — 0.2 sub-pixels. 

A 1.2 Marginalization 

We have experimented with marginalization over the amplitude 
components of the models. Since the model is linear in these pa- 
rameters their likelihood through a given slice in the other parame- 
ters is Gaussian and we can analytically marginalize over them. We 
have found this to be slightly problematic in general - the marginal- 
ization is not particularly numerically stable and since we have two 
components we need to perform the PSF convolution twice if we 
sum them after marginalization. These two effects meant it was in 
fact slower to use this trick, provided that we have a reasonable 
initial estimate of the ellipticity, particularly as our gradient-based 
minimizer finds the linear parameters very simple to optimize. 

If we move the modelling into Fourier space we could further 
marginalize over the centroid parameters xq and yo of the model - 
combining these dimensionality reductions could well yield a sig- 
nificant speed up. The downside of this is decreased flexibility in 
the modelling. We defer further discussion of this issue. 

A1.3 Point spread function (PSF) 

IM3SHAPE has the general facility to use an arbitrary image for 
PSF, but by default we read a catalog of either Moffat or Airy func- 
tion parameters and use that to build a PSF image. It is built at the 
same upsampled resolution as the image itself, and transformed and 
stored in Fourier space for faster convolution with trial images. 

When thes e images are made the co nvolution is done using the 
FFTW package iFrigo & Johnson] d2005h . 

A1.4 Fourier Sampling Kernel 

Our model for the pixels in an image is really a constant value 
across the surface of the pixel. But as far as a Discrete Fourier 
transform is concerned a 2D array represents a grid of 8— function 



samples from the image. To account for this difference our Fourier- 
space model of the image needs to include convolution with a 
square top-hat function the size of a pixel - this is true even when 
simulating at the same resolution as the data. Since we are also con- 
volving with a PSF we can fold together these convolutions into a 
single Fourier-space product. The Fourier transform of a 2D top-hat 
is given by: 

where the image dimensions are (n x X n y ) and the upsampling 
factor is u. 

A2 Resolution 

The model image in IM3SHAPEis made at a resolution several times 
higher than the data resolution, and the convolution with the PSF 
also happens at this high resolution. Only then is the model down- 
sampled to compare to the data. This procedure ensures that small 
and sharp features in the galaxy can be captured. 

We find, though, that the central regions of galaxy bulges with 
the standard de Vaucouleurs profile are so sharp that an even higher 
resolution is needed to accurately model them. Rather than generate 
the whole image at this extreme resolution we only use it for the 
central pixels where it is relevant. 

There are therefore three different resolution parameters to 
consider: the upsampling (the ratio of the model to image pixel size 
in the outer image region), the size of the central higher resolution 
region, and the resolution increase in this central region (compared 
to the outer region). The values that we used for these parameters 
are discussed in more detail in section |4~T1 



APPENDIX B: NUMERICAL OPTIMIZATION 

The optimization problem of finding the best model parameters 
can be tackled with a choice of algorithms; we have found that 
the best performing is the Levenberg-Marquadt (LevMar) method. 
This method uses more information than methods (such as Pow- 
ell's method or Nelder-Meald) which operate only on the total like- 
lihood as a function of the input parameters: it uses the complete 
mismatch map, (i.e. the \ 2 per pixel) for each parameter set con- 
sidered. 

The LevMar method uses gradients (analytic or numerical) to 
optimize; it is therefore not possible to signal a failure of a par- 
ticular parameter set (which can occur when parameters are near 
their limits) using NaNs or other flag values. While we do try to 
put in place a parameterization that prevents this from happening 
(see below), we also signal the failure by setting the \ 2 for the fail- 
ing model to 2 per pixel, a value large enough to be rejected but 
hopefully not hugely problematic. 

Bl Termination criteria 

The LevMar method depends on a number of manually set thresh- 
olds that determine when the numerical optimization should termi- 
nate. In particular, it stops and reports the parameters of the best 
model fit when at least one of the following criteria is met: 

• One of the partial derivatives of the objective function (x 2 ) 
with respect to the model parameters falls below a threshold si, 
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i.e. max < ex, where p denotes one of the model param- 

eters. Note that this termination criterion can be sensitive to the 
changes of image pixel scaling and changes in the signal-to-noise 
ratio (SNR). 

• The relative change of the norm of the estimated parameter 
vector is smaller than e 2 , i.e. V \p % ~ 1 -p'| 2 < e\ • b 4_1 | 2 . 
where i is the current iteration index. This termination criterion can 
be sensitive to individual parameter scaling. 

• The value of \ 2 falls below £3; again, this termination crite- 
rion is sensitive to image scaling and the S/N ratio. 

• The relative change of \ 2 fells below £4, i.e. |(x 2 ) 1-1 — 
(X 2 y\/(x 2 )'^ 1 < U./n, where (x 2 )' is the objective value at the 
current iteration and n is the total number of pixels. This termina- 
tion criterion was implemented by the authors and is not included 
in the LevMar implementation we are using. It can be sensitive to 
the step size that LevMar is taking and thus will be affected by the 
levmar_init_mu parameter - the initial damping of the normal 
equations solved in the Levenberg - Marquadt algorithm. 

• The maximum number of iterations is reached (the minimiza- 
tion has failed). 




Figure CI. Actions in the IM3SHAPEcodemn to prepare for the analysis of 
a single galaxy. Inputs, which can be collectively specified for a collection 
of galaxies, are upper red parallelograms, stored data are lower green par- 
allelograms, and calculations are yellow rectangles. The preparation time is 
negligible compared to the time of the optimization itself. 



APPENDIX C: ANALYSIS OVERVIEW 

Before analyzing a collection of galaxy images we: 

(i) Read a collection of options specifying the model to be used, pa- 
rameter ranges, and various options from a parameter file 

(ii) Build a model object encapsulating the model we will fit to the 
galaxies, including its likelihood and prior functions 

Before analyzing each individual image we then: 

(i) Generate the point-spread image based on PSF parameters; com- 
pute its FT. 

(ii) Generate the sine function for the specified upsampling (the FT of 
a top-hat). 

(iii) Multiply these two Fourier space quantities together to generate 
the Fourier kernel. 

(iv) Generate a weight map based on the image noise levels. 

(v) Compute a starting estimate for the galaxy properties from de- 
faults and weighted quadmpole moments. 

We then pass the model likelihood function to a minimizer 
object which runs the chosen algorithm. At each evaluation of the 
posterior, for our standard model, we: 

(i) Check if the parameters are outside the desired ranges and give up 
early if not. 

(ii) Build high-resolution (up-sampled) real-space models of each 
component of the galaxy (e.g. bulge and disc). 

(iii) Replace the central few pixels of each model with summed up 
pixels from even higher resolution, since this is the region where 
the image changes quickly. 

(iv) Sum and convolve the components with the precomputed kernel. 

(v) Down-sample the image to the data resolution. 

(vi) Compute the likelihood from the weight map and model image. 

(vii) Save the model image for use in the optimizer and residual analy- 
sis. 




Figure C2. The calculation performed whenever the optimization code pro- 
vides a set of parameters for which to calculate the image and likelihood. 
The principal parameterization used in this work and the normalization ap- 
plied to it are described in section [A] The "overwrite center" process is 
described in more detail in section ET] 



APPENDIX D: DETAILS ON THE MAXIMUM 
LIKELIHOOD ESTIMATION 

The maximization problem outlined in section [B] is based on the 
following model: 



y = D(K*x(p)) + n, 



(Dl) 



where y denotes the observed galaxy image, K = f *TI, the convo- 
lution of the PSF f and a square top-hat pixel window function II 
(see section lA1.4t . x(p) the (higher resolution) galaxy model with 



14 J. Zuntz et al. 



parameter vector p, and D(...) the downsampling operator, which 
averages into lower resolution pixels. As detailed in section [3] we 
assume a two-component Sersic galaxy model with the parameters 
p = (xq, 2/0, ei, e2,r d , A b , A d ) (see Table[2]for a detailed param- 
eter description) 

x (p) = ^2 A k exp|-R(a;o,J/o,ei,e2,r- d ) 1/2nfc | , (D2) 





Figure Dl. Components of Jacobian matrix with respect to its model pa- 
rameters (from left to right): xq, yo, e\, e 2 , r d , A b , A d . 



where = 1 for the disc component and 4 for the bulges, and R 
denotes the sheared squared distance 



Ri = d]Cdi 



(D3) 



with the centered coordinate vector d, = (xi — xo, y% — yo) and 
the covariance matrix 



C = 



1 



rd(l-|e|) 



1 + |e| — 2ei -2e 2 
-2e 2 1 + |e] + 2ei 



(D4) 



For additive Gaussian noise, i.e. n oc J\f(0, S) maximizing the 
likelihood P(p|y, K, x) is equivalent to minimizing 



(D5) 



with respect to the model parameters p, where r denotes the resid- 
ual image between the observed galaxy image y and the downsam- 
pled, PSF convolved galaxy model image x(p), i.e. 



r(p) 



[y-D(K*x(p))] 



(D6) 



For simplicity we assume n to be uncorrected white noise and so 
consider only diagonal terms of the covariance matrix, i.e. a — 
diag S; we neglect any correlation in noise of adjacent pixels. 

For minimizing dD5t , IM3SHAPE uses the free software 
package LEVMAbI 13 L a C/C++ implementation of the Levenberg- 
Marquardt (LM) algorithm, a standard techniqu e for solving non- 
linear least-squares problems dPress et al.|[l992l) . The LM method 
is based on a linear approximation of r in the neighbourhood of p, 



r(p + Ap) « r(p) + J(p)Ap + 0(]|Ap|| 



(D7) 



where Ap is a small perturbation of the parameter vector p and J 
denotes the Jacobian matrix of r: 



J «(P) = JT- 



(D8) 



The Jacobian dD8t can be numerically approximated, but to speed 
convergence of the minimization procedure it is generally recom- 
mended to compute it analytically. For our model (1DH the Jacobian 
reads with respect to the model parameters p £ {a;o , yo , ei , e 2 , } 



dr 

dp 



/2n k -l 



<9R 

dp 

(D9) 



where • denotes element-wise multiplication and the partial deriva- 
tives dH/dp are readily computed from tD3\ and (ID4t . With re- 
spect to the bulge and disc amplitudes A ke r b d \ one obtains 



^-Id(k.-p{-r^}) 



(D10) 



13 Available at |http : / /www . ics . forth . gr / -lourakis/levmar /| 



Note that for clarity and the ease of exposition we omit- 
ted the image upsampling within a central image region as de- 
scribed in section 14.11 in the derivation above, i.e. we consid- 
ered the special case where n_ central jpixeljap sampling is set 
to one. The generalization to n_central_pixel_up sampling > 1 
is straightforward since it only involves an additional summation 
over all n_centraLpixeLupsampling 2 sub-pixels to yield the flux 
[x(p)]jj for each pixel (i, j) within the central image region. 



APPENDIX E: NOISE BIAS CALIBRATION USING DEEP 
DATA 

Noise bias is the difference between the expectation of the 
maximum-likelihood result found by a model fitter like IM3SHAPE, 
which is a biased estimator, and the true underlying value, coming 
from the non-linear mapping between parameters and image. Its 
typical value is a few percent, and we need to remove it to reach 
our target accuracy levels. 

The size of the bias depends sensitively on the true parame- 
ters of the galaxy, and if these were known we could remove noise 
bias completely. However, we have access only to noisy estimates 
of these parameters and therefore our estimate of the noise bias is 
itself noisy, and biased. We can think of this as a bias-on-bias prob- 
lem, and we find that failing to account for it means we significantly 
miss our targets. 

We can get around this problem if there is a subset of our 
observational data that is deeper than the bulk of our sample, and 
therefore of greater signal-to-noise, to a degree sufficient that it has 
negligible noise bias. This is often the case in real surveys that seek 
to detect high redshift supernovae, and it is approximated in the 
GREAT08 challenge with LowNoise_Blind data set (although the 
galaxy types in this set do not match those in the main sample ex- 
actly). In this case we do not calibrate each galaxy individually, but 
instead find a mean bias for the population and apply it en masse. 

Biases m(9) and c(8) will afflict each of our galaxies, depend- 
ing on the true galaxy properties 9. We can calculate mean values 
of this bias m and c, and apply these evenly to all the galaxies, such 
that the mean galaxy ellipticity (which is our goal) is correct. 

The mean of the bias across the population is given by: 



rh= I m(8)p(9)A8. 



(El) 



We find the distribution p(8) using fits to the deep data, and m(9) 
using simulations - for each point 9 V in a grid in the parameter 
space we simulate many galaxies and determine the value m(8 p ). 
We can then interpolate between these values to do the integral, and 
finally apply the mean m to all the galaxy estimates. 

A similar process is used for the c bias, except that c varies 
directly with the PSF, so we use this information - we calculate the 
c assuming a fiducial PSF with ellipticity eg af aligned with the ei 
direction. The applied value to apply to each galaxy is then: 



Cl + IC2 



psf 



»psf 



(E2) 
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Table Fl. Quality factor Q scores for GREAT08 from this work and the 
original challenge entrants. 



Name 


1 


2 


3 


4 


Branch 
5 


6 


7 


8 


9 


13 


817.2 


1342.7 


450.8 


725.1 


722.0 


4702.5 


84.6 


382.7 


99.8 


I3-U 


128.7 


183.5 


67.2 


78.1 


115.4 


1417.1 


25.3 


2613.0 


10.6 


HB 


256.4 


654.4 


241.9 


184.2 


158.9 


225.3 


176.0 


272.1 


131.7 


AL 


764.9 


595.1 


348.4 


121.9 


404.0 


282.8 


22.4 


529.5 


512.1 


TK 


224.7 


569.6 


737.4 


161.9 


238.5 


63.6 


54.3 


142.9 


59.9 


CH 


174.2 


324.7 


157.4 


83.3 


212.4 


173.9 


13.1 


230.8 


18.6 


PG 


21.7 


26.5 


26.3 


37.8 


44.2 


28.5 


50.6 


23.1 


82.4 


MV 


246.2 


155.3 


95.1 


126.3 


224.4 


488.7 


4.3 


398.8 


23.8 


KK 


222.2 


554.8 


448.6 


89.8 


230.4 


100.0 


3.1 


163.4 


33.7 


HHS3 


26.9 


33.0 


30.7 


10.2 


20.3 


30.2 


30.3 


20.5 


25.2 


SB 


161.4 


93.4 


226.3 


114.5 


155.2 


19.8 


4.9 


26.0 


8.6 


HHS2 


23.0 


26.8 


26.1 


9.9 


18.3 


23.4 


25.4 


19.4 


23.2 


HHS1 


67.2 


109.0 


22.3 


43.0 


60.9 


232.7 


8.8 


1516.6 


2.1 


MJ 


10.6 


11.5 


10.4 


23.0 


4.0 


8.3 


13.9 


10.4 


16.4 


USQM 


13.3 


15.8 


11.4 


9.4 


23.4 


8.4 


0.8 


1.1 


0.2 



APPENDIX F: RESULTS TABLES 

Tables [FT1F41 give the GREAT08 Q, m, ci, and c 2 values for 
IM3SHAPE with (13) and without (I3-U) noise bias calibration ap- 
plied. Also shown are the score s for the other entra nts at the time 
of the challenge, as described in iBridle et al.N2010t) . 
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Table F2. Multiplicative bias m scaled by 10~ 3 on GREAT08 from this work and the original challenge entrants. 



Name 


1 


2 


3 


4 


Branch 
5 


6 


7 


8 


9 


13 


15.3 


3.3 


-3.3 


12.9 


13.9 


3.4 


-35.8 


-18.3 


91.9 


I3-U 


37.6 


25.4 


49.8 


35.1 


36.2 


8.9 


49.0 


3.3 


1 16.0 


HB 


28.9 


12.8 


25.6 


23.3 


24.9 


28.0 


11.3 


25.6 


37.8 


AL 


9.9 


-3.1 


7.8 


4.9 


6.5 


23.3 


-87.5 


-3.4 


8.4 


TK 


26.6 


14.3 


-5.6 


20.9 


21.1 


59.6 


-57.0 


37.0 


-3.4 


CH 


17.8 


7.9 


11.0 


11.7 


16.7 


9.5 


-35.1 


24.6 


0.1 


PC 


70.3 


58.1 


58.9 


70.5 


67.0 


38.3 


-38.6 


61.4 


42.1 


MV 


-27.4 


-35.8 


-43.2 


-36.6 


-30.0 


19.6 


-188.4 


-19.3 


-84.4 


KK 


25.7 


8.8 


5.0 


16.8 


21.6 


45.4 


-254.4 


27.6 


-61.0 


HHS3 


74.1 


64.1 


67.0 


70.2 


69.8 


66.2 


63.8 


87.6 


74.5 


SB 


-36.9 


-48.1 


-22.6 


-41.6 


-35.6 


-101.2 


-167.4 


-88.4 


-7.3 


HHS2 


83.1 


70.7 


72.7 


75.5 


79.2 


79.1 


75.8 


90.4 


82.5 


HHS1 


54.5 


39.1 


89.6 


43.7 


46.6 


-28.1 


150.7 


-8.9 


306.0 


MJ 


61.1 


47.1 


58.0 


54.9 


53.3 


-36.5 


22.5 


3.1 


26.6 


USQM 


-52.5 


-60.0 


-62.7 


-51.3 


-49.2 


-118.6 


420.8 


-393.5 


882.2 



Table F3. Additive bias c\ scaled by 10 5 on GREAT08 from this work and the original challenge entrants. 



Branch 

Name 123456789 



13 


-22.0 


-31.8 


63.8 


-33.9 


12.8 


-17.8 


-52.5 


44.6 


-254.3 


I3-U 


-67.0 


-77.0 


-94.6 


-124.0 


29.4 


-29.0 


-235.0 


1.1 


-304.0 


HB 


-4.4 


-17.6 


-16.8 


-66.3 


58.2 


-45.2 


59.4 


-26.1 


10.7 


AL 


-38.4 


-51.5 


-59.3 


-118.6 


34.8 


-54.5 


-79.9 


-53.3 


-51.3 


TK 


-26.3 


-28.1 


20.5 


-63.5 


9.5 


-50.7 


13.4 


-44.4 


154.6 


CH 


82.3 


61.7 


84.3 


143.2 


-22.0 


93.1 


-9.9 


23.7 


262.5 


PG 


188.4 


177.9 


182.1 


-41.9 


-87.5 


220.0 


-66.8 


186.5 


34.1 


MV 


8.5 


-0.9 


-37.5 


-31.6 


23.1 


-6.7 


-283.5 


38.0 


-79.3 


KK 


-49.2 


-44.6 


-47.1 


-132.2 


37.6 


-57.8 


-119.5 


-62.8 


-129.9 


HHS3 


-174.2 


-165.6 


-160.1 


-387.2 


75.8 


-180.0 


-168.8 


-191.8 


-196.5 


SB 


9.5 


16.6 


-54.5 


18.7 


16.0 


135.4 


-311.2 


110.3 


-459.1 


HHS2 


-177.6 


-185.8 


-176.0 


-385.1 


86.0 


-196.5 


-178.0 


-195.3 


-189.0 


HHS1 


-71.2 


-74.8 


-138.6 


-169.4 


56.4 


48.4 


-225.2 


1.2 


-431.7 


MJ 


-245.8 


-249.7 


-219.1 


69.1 


-696.9 


-345.8 


-92.1 


-288.8 


285.7 


USQM 


-138.0 


-111.2 


-114.0 


-295.7 


70.3 


-118.8 


-149.5 


-161.5 


-20.2 



Table F4. Additive bias C2 scaled by 10 5 on GREAT08 from this work and the original challenge entrants. 



Name 


1 


2 


3 


4 


Branch 
5 


6 


7 


8 


9 


13 


2.1 


-10.6 


-17.0 


9.1 


3.4 


-0.6 


9.4 


-24.4 


61.2 


I3-U 


18.5 


5.6 


41.7 


42.1 


47.9 


3.5 


75.8 


-8.5 


78.9 


HB 


15.3 


11.6 


31.2 


43.2 


28.1 


19.0 


31.6 


15.3 


-1.6 


AL 


19.6 


11.0 


39.7 


49.6 


49.3 


15.8 


-14.1 


10.7 


-1.0 


TK 


44.6 


30.8 


40.0 


69.8 


61.5 


12.5 


28.1 


24.5 


-92.7 


CH 


-16.2 


-37.5 


6.1 


-39.4 


-72.4 


-33.0 


53.4 


-37.7 


-24.0 


PG 


21.2 


14.7 


39.0 


81.9 


-9.0 


11.0 


79.9 


5.0 


40.3 


MV 


14.4 


-14.0 


26.9 


27.7 


-4.5 


-11.1 


82.0 


-1.5 


1.9 


KK 


18.7 


1.7 


41.8 


60.1 


39.4 


9.6 


10.2 


19.6 


-24.4 


HHS3 


54.9 


63.1 


83.7 


148.3 


193.6 


57.3 


70.9 


69.6 


56.5 


SB 


-15.1 


-27.0 


22.3 


-32.8 


-45.2 


-53.7 


87.5 


-43.4 


126.8 


HHS2 


69.4 


65.4 


84.5 


150.6 


185.9 


67.3 


58.1 


69.3 


48.3 


HHS1 


27.1 


18.2 


79.0 


72.4 


77.0 


-20.1 


49.2 


-4.7 


145.2 


MJ 


322.9 


316.9 


349.5 


218.3 


127.9 


313.0 


353.2 


329.2 


164.7 


USQM 


-114.4 


-54.5 


-182.7 


-83.3 


-33.3 


22.0 


-977.9 


187.4 


-1615.3 



